##### for Figure 4 e-j


rm(list=ls())

library(Seurat)
library(dplyr)
library(tidyverse)
library(patchwork)
library(devtools)
library(harmony)
library(clustree)
library(presto)
library(fgsea)


setwd("/Users/gin/fsdownload/mspleen_sc/datadir")

##### data format ##### -----
dir_name <- list.files()
names(dir_name) = c('stromal1',"stromal2")  

scRNAlist <- list()
for(i in 1:length(dir_name)){
  counts <- Read10X(data.dir = dir_name[i])
  scRNAlist[[i]] <- CreateSeuratObject(counts, min.features =200)
}


for(i in 1:length(scRNAlist)){
  sc <- scRNAlist[[i]]
  sc[["mt_percent"]] <- PercentageFeatureSet(sc, pattern = "^mt-")
  HB_genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
  HB_m <- match(HB_genes, rownames(sc@assays$RNA))
  HB_genes <- rownames(sc@assays$RNA)[HB_m] 
  HB_genes <- HB_genes[!is.na(HB_genes)] 
  sc[["HB_percent"]] <- PercentageFeatureSet(sc, features=HB_genes) 
  scRNAlist[[i]] <- sc
  rm(sc)
}

scRNAlist <- lapply(X = scRNAlist, FUN = function(x){
  x <- subset(x, 
              subset = nFeature_RNA > 500 & nFeature_RNA < 6000 & 
                mt_percent < 2 & 
                HB_percent < 3 & 
                nCount_RNA < quantile(nCount_RNA,0.97) & 
                nCount_RNA > 1000)})

scRNAlist <- merge(scRNAlist[[1]],y=c(scRNAlist[[2]]))


scRNAlist <- NormalizeData(scRNAlist) %>% 
  FindVariableFeatures(selection.method = "vst",nfeatures = 3000) %>% 
  ScaleData() %>% 
  RunPCA(npcs = 30, verbose = T)


scRNA_harmony <- RunHarmony(scRNAlist, group.by.vars = "orig.ident",project.dim = F)

scRNA_harmony <- FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:30)
scRNA_harmony <- FindClusters(object = scRNA_harmony,resolution = c(seq(0,1,.1)))
clustree(scRNA_harmony)
scRNA_harmony <- FindClusters(scRNA_harmony, resolution = 0.1)

scRNA_harmony <- RunTSNE(scRNA_harmony, reduction = "harmony", dims = 1:20)
scRNA_harmony <- RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:20)

##### singleR annotation ##### -----
library(SingleR)
library(celldex)

sce_for_SingleR <- GetAssayData(scRNA_harmony,slot = "counts")
clusters=scRNA_harmony@meta.data$seurat_clusters
mouseImmu <- ImmGenData()
pred.mouseImmu <- SingleR(test = sce_for_SingleR, ref = mouseImmu, labels = mouseImmu$label.main,
                          clusters = clusters, 
                          assay.type.test = "logcounts", assay.type.ref = "logcounts")

mouseRNA <- MouseRNAseqData()
pred.mouseRNA <- SingleR(test = sce_for_SingleR, ref = mouseRNA, labels = mouseRNA$label.fine ,
                         clusters = clusters, 
                         assay.type.test = "logcounts", assay.type.ref = "logcounts")

cellType=data.frame(ClusterID=levels(scRNA_harmony@meta.data$seurat_clusters),
                    mouseImmu=pred.mouseImmu$labels,
                    mouseRNA=pred.mouseRNA$labels)
cellType

scRNA_harmony@meta.data$cellType = "NA"
cellType$mouseRNA[2] = "P210+ Granulocytes"
cellType$mouseRNA[7] = "Fibroblasts"
cellType$mouseRNA[12] = "Fibroblasts"
for (i in 1:nrow(cellType)){
  scRNA_harmony@meta.data[which(scRNA_harmony@meta.data$seurat_clusters == cellType$ClusterID[i]),"cellType"] <- cellType$mouseRNA[i]
}
scRNA_harmony@meta.data$cellid = 1:ncol(scRNA_harmony)

DimPlot(scRNA_harmony,reduction = "tsne",group.by = "cellType",label=T)

##### fibro cluster ##### -----
fibro = subset(scRNA_harmony,cellType=="Fibroblasts")

fibro <- NormalizeData(fibro, normalization.method = "LogNormalize") 
fibro <- FindVariableFeatures(fibro, selection.method = 'vst', nfeatures = 2000)
fibro <- ScaleData(fibro)
fibro <- RunPCA(fibro, features = VariableFeatures(object = fibro)) 

fibro <- FindNeighbors(fibro, dims = 1:20)
fibro <- FindClusters(fibro, resolution = 0.1 )
# Look at cluster IDs of the first 5 cells

fibro <- RunUMAP(fibro, dims = 1:20)
fibro <- RunTSNE(fibro, dims = 1:6)
DimPlot(fibro, reduction = 'umap')
DimPlot(fibro, reduction = 'tsne',pt.size = 2)
levels(fibro)

FeaturePlot(fibro,"Grem1",reduction = "tsne",order=T,pt.size=2)
FeaturePlot(fibro,"Grem1",reduction = "tsne",order=T,pt.size=2,min.cutoff = 0)
VlnPlot(fibro,"Grem1")

seu_cluster = c("0","1","2","3")
fb_name = c("Fibroblast1","Fibroblast2","Fibroblast3","Fibroblast4")
fb_type = as.data.frame(cbind(seu_cluster,fb_name))

fibro@meta.data$fb_type = "NA"

for (i in 1:nrow(fb_type)){
  fibro@meta.data[which(fibro@meta.data$seurat_clusters == fb_type$seu_cluster[i]),"fb_type"] <- fb_type$fb_name[i]
}

fibro@meta.data$cellType = fibro@meta.data$fb_type


for (i in 1:ncol(fibro)){
  scRNA_harmony@meta.data[which(scRNA_harmony@meta.data$cellid == fibro@meta.data$cellid[i]),"cellType"] = fibro@meta.data$cellType[i]
}


DimPlot(scRNA_harmony,reduction = "tsne",group.by = "cellType",label = T)
DimPlot(fibro, reduction = 'tsne',pt.size = 2,group.by = "cellType",label = T)

##### p210 cluster ##### -----
p210 = subset(scRNA_harmony,cellType=="P210+ Granulocytes")

p210 <- NormalizeData(p210, normalization.method = "LogNormalize") 
p210 <- FindVariableFeatures(p210, selection.method = 'vst', nfeatures = 2000)
p210 <- ScaleData(p210)
p210 <- RunPCA(p210, features = VariableFeatures(object = p210)) 

p210 <- FindNeighbors(p210, dims = 1:20)
p210 <- FindClusters(p210, resolution = 0.1 )
# Look at cluster IDs of the first 5 cells

p210 <- RunUMAP(p210, dims = 1:20)
p210 <- RunTSNE(p210, dims = 1:10)
DimPlot(p210, reduction = 'umap')
DimPlot(p210, reduction = 'tsne',pt.size = 2)


seu_cluster = c("0","1","2","3")

marker = FindAllMarkers(p210,only.pos = T,logfc.threshold = 0.25,min.pct=0.1)
VlnPlot(p210,"P210")

seu_cluster = c("0","1","2","3")
p210_name = c("P210+_1","P210+_2","P210+_3","P210+_4")
p210_type = as.data.frame(cbind(seu_cluster,p210_name))

p210@meta.data$p210_type = "NA"

for (i in 1:nrow(p210_type)){
  p210@meta.data[which(p210@meta.data$seurat_clusters == p210_type$seu_cluster[i]),"p210_type"] <- p210_type$p210_name[i]
}

p210@meta.data$cellType = p210@meta.data$p210_type

#反向mapping
for (i in 1:ncol(p210)){
  scRNA_harmony@meta.data[which(scRNA_harmony@meta.data$cellid == p210@meta.data$cellid[i]),"cellType"] = p210@meta.data$cellType[i]
}
DimPlot(scRNA_harmony,reduction = "tsne",label = T,group.by = "cellType",pt.size = 1)

##### p210 monocle ##### -----

setwd("/Users/gin/fsdownload/mspleen_sc/workdir/")

library(monocle)

umi <- as(as.matrix(p210@assays$RNA@counts), 'sparseMatrix')

pData <- p210@meta.data
pData$celltype <- p210@active.ident

fData <- data.frame(
  gene_short_name = row.names(p210),
  row.names = row.names(p210)
)

pd <- new('AnnotatedDataFrame', data=pData)
fd <- new('AnnotatedDataFrame', data=fData)
cds <- newCellDataSet(
  umi,
  phenoData = pd,
  featureData = fd,
  lowerDetectionLimit = 0.1,
  expressionFamily = negbinomial.size()
)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

expressed_genes <- row.names(
  subset(fData(cds),)
)
diff_test_res <- differentialGeneTest(
  cds[expressed_genes,],
  fullModelFormulaStr = "~cellType")

ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))

cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)


library(Seurat)
express_genes <- VariableFeatures(p210)
cds <- setOrderingFilter(cds, express_genes)
plot_ordering_genes(cds)

deg.cluster <- FindAllMarkers(p210)
express_genes <- subset(deg.cluster, p_val_adj < 0.05)$gene
cds <- setOrderingFilter(cds, express_genes)
plot_ordering_genes(cds)

disp_table <- dispersionTable(cds)
disp.genes <- subset(disp_table, mean_expression >= 0.1 &
                       dispersion_empirical >= 1*dispersion_fit)$gene_id

cds <- setOrderingFilter(cds, disp.genes)
plot_ordering_genes(cds)
cds <- reduceDimension(cds, max_components = 2,
                       method = 'DDRTree')

cds <- orderCells(cds)
plot_cell_trajectory(cds)


cds$tmp = cds@assayData$exprs
plot_cell_trajectory(cds, color_by = "Pseudotime")
plot_cell_trajectory(cds, color_by = "cellType")
cds <- orderCells(cds, root_state = 3)
plot_cell_trajectory(cds)


##### Grem1+ marker gene ##### -----
DimPlot(scRNA_harmony,reduction = "tsne",group.by = "cellType",pt.size = 1,label = T)
FeaturePlot(scRNA_harmony,"Grem1",reduction = "tsne",order = T,pt.size = 1)
DimPlot(fibro,reduction = "tsne",group.by = "cellType",pt.size = 2,label = T)
VlnPlot(fibro,"Grem1",pt.size=0,group.by = "cellType")
FeaturePlot(fibro,"Grem1",reduction = "tsne",order = T,pt.size = 2)

marker = FindMarkers(scRNA_harmony,ident.1 = "Fibroblast4",
                     group.by = "cellType",only.pos = T,
                     logfc.threshold = 0.25,min.pct = 0.2)
write.csv(marker,"allGrem1marker.csv")


##### heatmap ##### -----
library(scRNAtoolVis)
library(Seurat)

genelist = c("Ly6a","Tlx1","Wt1","Tcf21",#RP cluster
             "Gsn","Dpt","Thy1","Il33", #TRC cluser
             "Mcam","Rgs5","Cspg4","Cnn1",#mural cells express the classical pericyte markers
             "Cxcl13","Madcam1","Ch25h","Igfbp3","Grem1") #MRC cluster

jjDotPlot(object = fibro,
          gene = genelist,id = "cellType")

marker = FindAllMarkers(fibro,only.pos = T,logfc.threshold = 0.25,min.pct=0.1)

AverageHeatmap(object = fibro,
               gene = marker$gene,
               clusterAnnoName = F,
               showRowNames = F,
               markGenes = genelist,fontsize = 6,
               group.by = "cellType")

##### GSEA ##### -----

library(presto)

genes = wilcoxauc(fibro,"cellType")
gsea_genes<-genes %>% dplyr::filter(group=="Fibroblast4") %>% arrange(desc(auc)) %>% dplyr::select(feature,auc)
ranks = deframe(gsea_genes)


library(msigdbr)
m_df = msigdbr(species = "Mus musculus", category = "C5")
fgsea_sets = m_df %>% split(x = .$gene_symbol, f = .$gs_name)

fgsea_res = fgsea(fgsea_sets, stats = ranks, nperm = 1000)

fgsea_tidy = fgsea_res %>% as_tibble() %>% arrange(desc(NES))

fgsea_tidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% head()

write.csv(fgsea_tidy[,1:7],"GSEA_C5.csv")

m_df = msigdbr(species = "Mus musculus", category = "C2")
fgsea_sets = m_df %>% split(x = .$gene_symbol, f = .$gs_name)

fgsea_res = fgsea(fgsea_sets, stats = ranks, nperm = 1000)

fgsea_tidy = fgsea_res %>% as_tibble() %>% arrange(desc(NES))

fgsea_tidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% head()

write.csv(fgsea_tidy[,1:7],"GSEA_C2.csv")

m_df = msigdbr(species = "Mus musculus", category = "C3")
fgsea_sets = m_df %>% split(x = .$gene_symbol, f = .$gs_name)

fgsea_res = fgsea(fgsea_sets, stats = ranks, nperm = 1000)

fgsea_tidy = fgsea_res %>% as_tibble() %>% arrange(desc(NES))

fgsea_tidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% head()

write.csv(fgsea_tidy[,1:7],"GSEA_C3.csv")


##### bubble plot ##### -----
library(ggplot2)
bubble = read.table("selected GSEA.txt")
bubble = bubble[,c(1,2,5,7)]
bubble[,2] = -log10(bubble[,2])
#bubble[,4] = log10(bubble[,4])


labels=bubble[order(bubble$NES),"pathway"]
bubble$pathway = factor(bubble$pathway,levels=labels)

ggplot(bubble,aes(x=NES,y=pathway))+
  geom_point(aes(size=size,color=pval))+ 
  scale_colour_gradient(low="orange",high="red",)+
  scale_size_continuous(range = c(12,18)) + 
  xlim(c(2,18))+
  theme_bw()

DimPlot(scRNA_harmony,reduction = "tsne",group.by = "cellType",label=T)
DimPlot(p210,reduction = "tsne",group.by = "cellType",label=T)

##### SingleCellSignalR ##### -----
#BiocManager::install("SingleCellSignalR")
library(SingleCellSignalR)
DimPlot(scRNA_harmony,group.by = "cellType",reduction = "tsne")
p4 = subset(scRNA_harmony,subset = cellType == "P210+_4")
f4 = subset(scRNA_harmony,subset = cellType == "Fibroblast4")

data1 = as.data.frame(p4[["RNA"]]@counts)
data2 = as.data.frame(f4[["RNA"]]@counts)
data = cbind(data1,data2)

genes = rownames(data)
cluster = cluster=c(rep(1,length(colnames(p4))),rep(2,length(colnames(f4))))
signal = cell_signaling(
  data,
  genes,
  cluster,
  #  int.type = c("paracrine", "autocrine"),
  c.names = NULL,
  s.score = 0.5,
  logFC = log2(1.5),
  species = "mus musculus",
  tol = 1,
  write = TRUE,
  verbose = TRUE
)
#1-2 P210+ Fibroblast4
#2-1 Fibroblast4 P210+

visualize_interactions(
  signal,
  #  show.in = NULL,
  #  write.in = NULL,
  write.out = FALSE,
  method = "default",
  limit = 30
)
signal1 = signal[[1]]
signal1$f4 = rep(0,length(rownames(signal1)))
signal1$p4 = rep(0,length(rownames(signal1)))
i=1
while (i <= length(rownames(signal1))){
  gene = print(signal1$`cluster 1`[i],quote=F)
  pct = length(which((data1[gene,]>0)))/length(colnames(p4))
  signal1$f4[i] = pct
  i = i+1
}
i=1
while (i <= length(rownames(signal1))){
  gene = print(signal1$`cluster 2`[i],quote=F)
  pct = length(which((data2[gene,]>0)))/length(colnames(f4))
  signal1$p4[i] = pct
  i = i+1
}
signal1 = signal1[signal1$f4>0.15&signal1$p4>0.15,]

signal2 = signal[[2]]
signal2$f4 = rep(0,length(rownames(signal2)))
signal2$p4 = rep(0,length(rownames(signal2)))
i=1
while (i <= length(rownames(signal2))){
  gene = print(signal2$`cluster 2`[i],quote=F)
  pct = length(which((data2[gene,]>0)))/length(colnames(f4))
  signal2$f4[i] = pct
  i = i+1
}
i=1
while (i <= length(rownames(signal2))){
  gene = print(signal2$`cluster 1`[i],quote=F)
  pct = length(which((data1[gene,]>0)))/length(colnames(p4))
  signal2$p4[i] = pct
  i = i+1
}
signal2 = signal2[signal2$f4>0.15&signal2$p4>0.15,]
write.csv(signal1,"signal1.csv")
write.csv(signal2,"signal2.csv")
tmp = list(signal1,signal2)
names(tmp) = c("1-2","2-1")
visualize_interactions(tmp,show.in = "1-2")
visualize_interactions(tmp,show.in = "2-1",limit = 15)



##### cellchat ##### -----
library(CellChat)
#提取表达矩阵和细胞种类信息
data.input <- GetAssayData(scRNA_harmony, assay = "RNA", slot = "data")
identity <- subset(scRNA_harmony@meta.data, select = "cellType")
#创建cellchat对象
cellchat <- createCellChat(object = data.input, meta = identity,  group.by = "cellType") #注意group.by参数
#导入数据库
####可选CellChatDB.human, CellChatDB.mouse
CellChatDB <- CellChatDB.mouse
#可以查看下数据库的信息
showDatabaseCategory(CellChatDB)
########在CellChat中，我们还可以先择特定的信息描述细胞间的相互作用，可以理解为从特定的侧面来刻画细胞间相互作用，比用一个大的配体库又精细了许多。
##查看可以选择的侧面都有哪些
unique(CellChatDB$interaction$annotation)
# use Secreted Signaling for cell-cell communication analysis，选择从想要研究的方面来特定筛选
CellChatDB.use <- subsetDB(CellChatDB, search = "Cell-Cell Contact")
CellChatDB.use <- CellChatDB
#赋值进去
cellchat@DB <- CellChatDB.use # set the used database in the object

###对表达数据进行预处理
##将信号基因的表达数据进行子集化，以节省计算成本
cellchat <- subsetData(cellchat)

# 识别过表达基因
cellchat <- identifyOverExpressedGenes(cellchat)
# 识别配体-受体对
cellchat <- identifyOverExpressedInteractions(cellchat)
# 将配体、受体投射到PPI网络
cellchat <- projectData(cellchat, PPI.mouse) #PPI.mouse

##相互作用推断
## 1、计算通信概率推断细胞互作的通信网络
cellchat <- computeCommunProb(cellchat, raw.use = F)
###如果特定细胞群中只有少数细胞，则过滤掉细胞间的通信，要筛选在多个细胞中有的互作
cellchat <- filterCommunication(cellchat, min.cells = 3)

#提取推断出的细胞互作的通信网络数据框，我们提供了一个subsetCommunication 函数，
#可以方便地访问感兴趣的推断的细胞间通信。
##返回一个数据框，包含所有推断的配体/受体级别的细胞-细胞通信。设置slot.name = "netP"以访问信令路径级别的推断通信
df.net <- subsetCommunication(cellchat) #返回的结果就是计算结果
levels(cellchat@idents)            #查看细胞顺序

##2、在信号通路水平上推断细胞间的通讯，看细胞间整体效应
cellchat <- computeCommunProbPathway(cellchat)
##汇总通信概率来计算细胞间的聚合通信网络。
cellchat <- aggregateNet(cellchat)
##3、计算聚合细胞互作通信网络
groupSize <- as.numeric(table(cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")
#左图：外周各种颜色圆圈的大小表示细胞的数量，圈越大，细胞数越多。发出箭头的细胞表达配体，
#箭头指向的细胞表达受体。配体-受体对越多，线越粗。
#右图：互作的概率或者强度值（强度就是概率值相加）



##每个细胞如何跟别的细胞互作（互作的强度或概率图）——weight
mat <- cellchat@net$weight
#par布局
par(mfrow = c(1,1), xpd=TRUE)

for (i in 1:nrow(mat)) {
  pdf(paste0(print(i),".pdf"))
  mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
  mat2[i, ] <- mat[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat), title.name = rownames(mat)[i])
  dev.off()
}

##每个细胞如何跟别的细胞互作（number+of+interaction图）——count
mat <- cellchat@net$count
par(mfrow = c(1,1), xpd=TRUE)
for (i in 1:nrow(mat)) {
  pdf(paste0(print(i),".pdf"))
  mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
  mat2[i, ] <- mat[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat), title.name = rownames(mat)[i])
  dev.off()
}

####可视化每个信号通路，针对于特定细胞通路
##查看通路

levels(cellchat@idents)            #查看细胞顺序
vertex.receiver = c(1, 2)          #指定靶细胞的索引
cellchat@netP$pathways             #查看富集到的信号通路
pathways.show <- "IL6"            #指定需要展示的通路


##层次图
vertex.receiver = seq(1,4) # a numeric vector. 
test = netVisual_aggregate(cellchat, signaling = "IL6",  vertex.receiver = vertex.receiver,layout="hierarchy")
#在层次图中，实体圆和空心圆分别表示源和目标。圆的大小与每个细胞组的细胞数成比例。线越粗，互作信号越强。看的是权重
#左图中间的target是我们选定的靶细胞。右图是选中的靶细胞之外的另外一组放在中间看互作。
ggsave("umap_tsne_integrated.pdf",test)
##圈图
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling ="IL6", layout = "circle")

##和弦图
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling ="IL6", layout = "chord", vertex.size = groupSize)

##热图
par(mfrow=c(1,1))
pdf("heatmap_weight.pdf")
netVisual_heatmap(cellchat, color.heatmap = "Reds", measure = c("weight"))
dev.off()

pdf("heatmap_count.pdf")
netVisual_heatmap(cellchat, color.heatmap = "Reds", measure = c("count"))
dev.off()
##纵轴是发出信号的细胞，横轴是接收信号的细胞，热图颜色深浅代表信号强度。
##上侧和右侧的柱子是纵轴和横轴强度的累积



